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Numerical Solution of Transonic Full Stream 

"k 

Function Equations in Conservation Form 


Abstract 


The stream function equation in conservation form is solved 
iteratively based on the artificial compressibility method. The 
density is not a unique function of the mass flux. To avoid the 
ambiguity near the sonic line, the density is updated in terms of 
the velocity, which is obtained through a simple integration of a 
first order equation step-by-step in the flow field. Iteration 
algorithms and finite difference approximations are discussed, and 
numerical results of both conservative and nonconservative calculations 
are presented. 


This work is sponsored by NASA Ames Contract No. 


NASA-10252. 



Flow Research Note No. 178 
October 1979 


-ii- 


Table of Contents 

Page 

Abstract i 

Table of Contents ii 

1. Introduction 1 

2. Stream Function Formulation 2 

3. Artificial Viscosity 4 

4. Shock Jump Conditions 5 

5. Artificial Compressibility 6 

6. Iteration Algorithms 8 

7. Mixed Flow Calculations 10 

8. Relaxation Method 12 

9. Numerical Results 13 

Appendix: Transonic Small Disturbance Equation 14 

References 16 

Figures 



Flow Research Note No. 178 
October 1979 


- 1 - 


1. Introduction 

The stream function formulation is suitable for design problems and 
internal flow as well as rotational flow calculations. The stream lines 
and potential lines (in the case of irrotational flows) provide an 
orthogonal grid with a natural body fitted coordinate. 

There are, however, some difficulties in solving the stream 
function equation: (1) the density is not a unique function of the mass 

flux, and (2) stream function exits only for two-dimensional and axi- 
symmetric flows. Extension to three-dimensional flows in terms of more 
than one stream function is not simple to implement. Due to these 
stumbling blocks, the stream function has not been recently used. 

In 1944 Emmons (Reference 1) solved the stream function equation by 
hand relaxation where the shock was fitted. In 1968 Sells (Reference 2) 
calculated only subcritical flows using line relaxations. Here transonic 
flows are calculated where shocks are captured, based on the artificial 
compressibility method (Reference 3), and the density is updated in terms 
of the velocity without ambiguity. 
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2. Stream Function Formulation 

For irrotational flows, the governing equations are: 


(pu) x + (pv) y = 0 


(1) 


u - v =0 

y x 


( 2 ) 




2.2 2 r}y - 1 

M (u + v - 1) 

00 ' 


(3) 


A potential function exits and it is governed by: 


(p<J> ) + (pcf> ) =0 

x x y y 


(4) 


p = 


(M 2 a 2 ) T " 

oo ' 


1 = [l - X - ~- 1 M 2 (cJ ) 2 + <\) 2 
2 <» VT x r y 


- 1 ) 


Y - 1 


(5) 


In terms of natural coordinates s , n Equation (4) becomes: 


(1 - M )4> + <j> = 0 . 

ss nn 


( 6 ) 


Similarly, a stream function exists such that: 




= <f> = v 

p y 


K 

— ^ = (j) = u 

P X 


(7) 


(;!* , fn . ) 

' p ^n ’ p ^s ' 


( 8 ) 
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Hence 



or (1 - M 2 )^ + = 0 . (11) 

ss nn 

While the density (the speed of sound) is uniquely determined in terms of 
2 2 2 

the speed q = u + v , from Equation (3), Equation (10) indicates that 
there are two values of the density for certain mass flux (Sells, Reference 
2, obtained rational approximations for the subsonic and supersonic 
branches) . 

We notice also that Equation (11) is a nonlinear mixed-type equation 
similar to Equation (6). 
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3. Artificial Viscosity 

It is well known (Reference 3) that Equation (6) can be solved by 
adding an artificial viscosity term either explicitly or indirectly 
(through upwind differencing), namely, 


(1 - M )<f> + 4> = -£<(> 


ss 


nn 


sss 


( 12 ) 


The corresponding stream function equation is 


(1 - m 2 h 

ss 


+ = —Gib 

nn sss 


(13) 


Equation (13) is indeed equivalent to Equation (12) with the augmented 
artificial viscosity term. 

Therefore, upwind differencing of Equation (11) in exactly the same 
way that Equation (6) is solved will yield nonconservative calculations. 
In the first case, however, the flow is irrotational and due to the 
nonconservative differencing, a proper mass balance across the shock is 
not achieved and the error can be represented as a source distribution at 
the shock. On the other hand, the mass is conserved automatically in the 
stream function formulation and due to the nonconservative differencing, 
a vorticity distribution at the shock is produced. 
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4. Shock Jump Conditions 

Equation (1) admits a discontinuous solution such that 


t pu ! ■ (37)5 M ■ 0 <«) 

(af) s M + I"]) - 0 as) 


In terms of the velocity potential. Equations (14 and 15) read 


- (f l t* y ]i ■ 0 

M = 0 


(16) 


(17) 


and in terms of the stream function, they are: 


M = 0 


(18) 


(§)i^HT] = 0 


(19) 


Either the shock is fitted according to the above shock jump relations, or 
it is captured through a proper conservative finite differencing, such 
that in the limit of small mesh, the shock jump conditioning are 
satisfied. 
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5. Artificial Compressibility 

In Reference (3), the full potential equation in conservation form is 
solved. An artificial viscosity is added through modifying the density, 
namely, 


(p<J) ) + (ptj) ) =0 

x x y y 


where 


p - p - M p s i s 


P A 
s s 


- p A + - p 
q X x q y y 


( 20 ) 

( 21 ) 

( 22 ) 


Centered differences are used everywhere in Equation (20) . Only p g is 
evaluated using upwind differencing. 

Similarly, the artificial viscosity can be added to the stream 
function equation as follows: 


Hence, 


or 


y , u . 
u = — ■ *- + y — pA 
p p 'x x 


-ip 

X V . 

v = — - y — p A 

p p y y 


(tI + Cfl - ° 


U = ily/p 


(23) 

(24) 

(25) 

(26) 


v = - if* x /p 


(27) 
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(28) 
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6. Iteration Algorithms 

In the potential claculations, given the density from previous 
iterations, the velocity potential is calculated from Equation (20). Then 
the density is updated using Equation (5) . The iteration loop is shown in 
Sketch 1. 



Sketch 1 

A similar iteration loop is implemented for stream function calculations 
as shown in Sketch 2. The algorithm converges only for subcritical flows. 



Sketch 2 

The iteration algorithm described in Sketch 2 is based on the following 
procedure: 
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new 


1 - 


ILjl 


' P old 


Y - 1 


(29) 


2 2 

Notice that p decreases if the mass flux (\L» + ip ) increases. This 

new T x y 

is true only for subsonic flows. Equation (10) can be rewritten in the 
form: 


F(p) = p Y+1 - ( 1 + M* )p 2 + MW + = 0 


(30) 


A Newton iteration on the single nonlinear equation in p 
(Equation (30)) is given by: 


J(p old ) ' <p „ew - Pold' ‘ - F(P old ) 


where 


J(p) = (y + l)p Y - 2p(l + M^) 


(31) 


Notice that the Jacobian is negative or positive depending on whether 

& A 

P q ^ is larger than or smaller than p (p at a = a ) . For pure 

subsonic or supersonic flows, Equation (31) may be used. But for mixed 
flows, perturbing the sonic conditions yields two possible solutions, and 
there is no obvious way to exclude anyone of them based on local 
considerations, unless the sonic line is fitted. 
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7 . Mixed Flow Calculations 

To avoid the ambiguity in the density updating, a different approach 
is adopted to calculate the speed given the stream function without 
referring to the old density. Starting with Equation (2): 



(32) 


(33) 


(34) 


Knowing u on a curve crossing the characteristics of Equation (33), u 
can be calculated throughout the flow field step by step. The conditions 
along such a curve could be either pure subsonic or pure supersonic, where 
u can be obtained using a Newton iteration. 

Once u is obtained, V is evaluated from Equation (35) 


v = - 


l|< 

_> 


( 35 ) 


The iteration loop is shown in Sketch 3. 
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8. Relaxation Method 

Horizontal line over relaxation (with additional d> term. Reference 

xt 

3) of the partial differential equation for ip , given p , is found to be 
fast. To explain this, consider the transonic small disturbance system of 
equations : 


v = (u - K)u - Ew 
y xx 


u = v 
y x 


\u lC J ,t/' 


\> ri 

, CC-’ 


aw = - u + w 

y x 


(36) 


The characteristics of system (36) is given by 


n 


Tv 


c (‘. 


X - (u - k)X - e/a = 0 


(37) 


In the limit — 0 , the system (36) is parabolic with y = constant as 

a triple characteristic. Similarly, the streamlines are the characteristics 
of Equation (12). 

Moreover, the iterative matrix of the horizontal over relaxation 
procedure seems to have a dominant eigenvalue allowing for even more 
acceleration. 
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9. Numerical Results 

Transonic flows over NACA 0012 airfoil have been calculated based 
on the full potential equation as well as the stream function equation. In 
these calculations, linearized boundary conditions are used. The stream 
function equation is solved using the iteration procedure described in 
Sketch 3. Both conservative and nonconservative results are presented. 

The iteration history of both vertical and horizontal line relaxations are 
plotted. Results for different magnitudes of the artificial viscosity in 
the artificial density formulation are also shown. 

Finally, we would like to mention that rotational flows can be 
easily calculated by modifying Equation (2) where the righthand side 
becomes the vorticity U)(i|0 . The density formula is also slightly 

modified by e where sOJO is the entropy. Also, other iterative 

procedures (i.e., ADI, Fast Solver, MAD, etc.) are, of course, applicable. 
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Appendix 

Transonic Small Disturbance Equation 

Chin and Rizzetta (Reference 4) solved the small disturbance 
equation: 


(1 - M' )tp + 
xx 


i|> =0 

yy 


(A-l) 


where 


1 - M 2 = 1 - M 2 - (y 

00 


+ 1)M u 

00 


Following Emmons (Reference 1), the small perturbation approximation 
of Equation (23), page 19, 


is 


y 



(A-2) 


Equation (A-2) is consistent with the approximation 


hence 
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or 


4 > = - ^ 

xy 


xx 


(A-3) 


Chin's calculations are nonconservative, not only because of Equation (A-l) 
but also because of Equation (A-2) . 
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Figure 3-1. Iteration history of vertical line relaxation 
W = 1.8 , = 0.85 . 
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Figure 3-4. Iteration history of vertical line relaxation 
a) = 1.8 , M = 0.98 . 
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Figure 4-1. 


Iteration history of horizontal line relaxation 
w = 1.8 , M = 0.85 . 
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■3. Iteration history of horizontal line relaxation 
co = 1.8 , M = 0.95 . 
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Figure 4-4. 


Iteration history of horizontal line relaxation 
w = 1.8 , M = 0.98 . 
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Figure 5-2. u versus 


M = 0.85 , viscosity factor =1.0 
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Figure 5-3. 


versus 


M = 0.85 , viscosity factor = 0.5 
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Figure 5-4. u versus x , = 0.9 . 
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Figure 5-5. u versus x 
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Figure 5-6. u versus x , M ro = 0.98 . 
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